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We present an alternative approach to setting initial data in general relativity. We do not use 
a conformal decomposition, but instead express the 3-metric in terms of a given unit vector field 
and one unknown scalar field. In the case of axisymmetry, we have written a program to solve 
the resulting nonlinear elliptic equation. We have obtained solutions, both numerically and from a 
linearized analytic method, for a perturbation of Schwarzschild. 
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I. INTRODUCTION 

Interferometry antennas (e.g. LIGO Q]) all over the 
world are going online to measure very fine dynamical 
variations of spacetime curvature known as gravitational 
waves. It is hoped that this will serve as a further evi- 
dence concerning Einstein's theory of general relativity. 
Eventually it might lead to a very new kind of observa- 
tional astronomy because gravitational waves emitted by 
astronomical sources reach us more or less unattenuated 
and give us information when other astronomical sources 
block the electromagnetic spectrum. Our understand- 
ing of gravitational waveforms emitted by various events 
is limited in the strong field case, particularly with pro- 
cesses involving black holes. These are important because 
they are likely to be powerful and therefore detectable Q • 
So for some years there has been a considerable effort to 
solve the binary black hole collision problem numerically. 

In order to solve Einstein's equations numerically one 
needs initial data, i.e. the description of the state of the 
system on an initial Cauchy surface. The usual approach 
is to use the conformal method, in which it is assumed 
that the initial data is conformally related to a reference 
metric. Often the reference metric is flat because the 
problem is then much simpler, although other cases have 
been investigated (see below) . The literature on the sub- 
ject is quite extensive; we do not discuss it in any detail, 
but instead refer to a recent review Q. The motiva- 
tion for using conformally flat initial data as done by e.g. 
Misner Q and Bowen and York [I| , is mathematical con- 
venience rather than physics. As was pointed out in [|| 
conformally flat slices of the Kerr metric are not likely to 
exist. Even in the case of non-spinning black holes it was 
shown that astrophysically correct initial data cannot 
be represented this way. 

Conventional wisdom holds that the unphysical part 
of an initial data set will be radiated away from the sys- 
tem, within the order of the light-crossing time. While 
this may well be true, the hypothesis has not been well 
tested. In order to do so, initial data sets constructed by 
different methods are needed. Of course such methods, 
including the one presented here, will involve some form 
of mathematical convenience. Provided the methods are 



really different, the influence of the precise form of the 
intial data on the subsequent spacetime evolution, can 
be investigated. 

The construction of black hole initial data that is not 
conformally flat, has been accomplished in recent years. 
The conformal method with a non-flat reference metric 
in Kerr-Schild form has been used to construct black hole 
data 

HUGS; see also[n|. The Kerr-Schild ansatz which 
was introduced in [T^ was a promising ansatz that is very 
different to the conformal approach. However, only the 
case of a particular perturbation was solved i n Il2l . and 
subsequent developments have been limited [13j. The 
motivation for the present work was to find a variation 
of the ansatz of ^2] f° r which there is a prospect of being 
able to find a solution in a variety of circumstances. The 
ansatz presented here satisfies this condition. It leads 
to a (nonlinear) elliptic equation in one unknown, and 
we have written a program to solve the equation. The 
initial data equation has been solved for a perturbation 
of Schwarzschild, using both the numerical program as 
well as an analytic method that linearizes the problem. 

Much use is made of computer algebra in this work, 
and we have used both Maple and Mathematica so as to 
provide a check on our calculations. 

We first present our ansatz and the constraint equa- 
tions (Sec.nj. Next, in Sec. IIIII we express initial data 
for the Schwarzschild geometry in terms of our ansatz. 
In Sec. IIVI we use our ansatz to find initial data for a 
perturbation of Schwarzschild, solving the problem both 
numerically and by linearization. Sec. IVI summarizes our 
results and also discusses possible further work. 



II. MODIFICATION OF THE KERR-SCHILD 
ANSATZ 

We attempted to modify the Kerr-Schild ansatz so that 
the constraint equations would reduce to a single elliptic 
equation - because then, under many circumstances, a 
solution to the problem can be expected to exist. The 
simplest way to do this is to require that the extrinsic 
curvature be zero 
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(2.1) 
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and the momentum constraints are automatically satis- 
fied. The Hamiltonian constraint reduces to 



= 0, 



(2.2) 



and we used computer algebra to investigate the nature 
of Eq. (|2.2I) for various forms of the 3-metric that involve 
a unit vector field and one scalar function in analogy with 
[l2| . The following ansatz seems to be acceptable 
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1 - 2V 



(2.3) 



with dij the Euclidean (flat) 3-metric, V a scalar field, 
and ki the normalized gradient of another scalar field 
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(2.4) 



It is useful to note a number of points about the ansatz 
defined by Eqs. (23), (23) and (|2~3)) : 

• Under the circumstances described later, the ansatz 
leads to a nonlinear elliptic equation for which a so- 
lution can be explicitly constructed. Whether the 
ansatz leads to solutions in more general circum- 
stances is, of course, unknown. 

• The assumptions (|2.1|) and l|2.4|l . were made to sim- 
plify the problem. To what extent these assump- 
tions can be relaxed, needs to be investigated. In 
order for the method to be astrophysically useful, 
some relaxation will be necessary so as to be able 
to allow for situations where there is momentum 
and spin. 



We have evaluated the Hamiltonian constraint in 
the case of axisymmetry, i.e. x l — (r, 9, ip), — 
diag(l, r 2 , r 2 sin 2 9), and with V and $ depending only 
on r and 9) . The reason for the restriction to axisymme- 
try was problem simplification, while at the same time 
including non-trivial situations such as two black holes 
and a perturbed Schwarzschild black hole. We used com- 
puter algebra to find 

c2oV rr + cnV r e + c^Vfie + c 5 (V r ) 2 + c 4 U r Ue 

+ c 3 (V e ) 2 +c 2 V e + c 1 V r + c V = 0. (2.5) 

The coefficients are functions of r, 9, V, $ and higher 
derivatives of Some of the coefficients are very long 
expressions, and so here we give only the coefficients of 
the principal part 
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The determinant of the principal part is 
1 , 8 



A = C20C02 
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r' 2 (1 -2V(r,9)) 



showing that the equation is elliptic provided V < 1/2. 



III. SCHWARZSCHILD GEOMETRY 

We now find an explicit representation of the 
Schwarzschild geometry that satisfies the ansatz ()2.3|l . 
This will be used for a number of purposes, including the 
setting of boundary data in more general situations, and 
being the zeroth order solution about which a perturbed 
Schwarzschild solution will be constructed. 

Consider the Schwarzschild 3-metric in isotropic coor- 
dinates (f, 9, ip) 

ds 2 = (l + (df 2 + f 2 d9 2 + f 2 sin 2 9dcp 2 ) (3.1) 

with the extrinsic curvature of the spacelikc initial slice 
Kij = 0. We can obtain the form (|2.3|l by requiring 



1 + — ) dr = dr 

2r . 



(3.2) 



which is easily integrated to give the coordinate transfor- 
mation 



2r m m 
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(3.3) 



where the integration constant has been chosen so that 
the event horizon is at r = f = m/2. Then in (r,9,tp) 
coordinates the metric satisfies ansatz (|2.3ll with 



* = ;■ nr) = V S {r) = \(l-{l + f_) y2 

On the event horizon at r = m/2 we have 
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(3.4) 
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and also we can make an asymptotic expansion to find 

2r 



lim V(t') 
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(3.6) 



As an additional check on our calculations, we substi- 
tuted Eq. (23) into Eq. l|23)l . obt aining = as ex- 
pected. 



IV. PERTURBATION OF THE 
SCHWARZSCHILD METRIC 

We now explore the case of a perturbed single 
Schwarzschild black hole with mass m = 1. We set 
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$ = 1/r and use the simple inner boundary condition 

v(r=\, =H + e p n (cos0) (4.1) 

to perturb the black hole where P n is the ro-th Legendre 
polynomial. In principle, the size of e is only limited by 
the fact that V < 1/2 to preserve ellipticity of Eq. 12. 51 1. 
so 



At the analytic level, the outer boundary condition is 
V(r,9) — > Vs(r) as r — > oo, but computationally we set 
V(r,G) to Vs(r) at r = roB which makes physical sense 
when roB is large. 

We solve the problem in two ways. First, we linearize, 
that is we ignore all terms in Eq. 1)2. 5fl of order e 2 , and we 
are able to write the solution as an infinite sum of eigen- 
functions. Second, we develop a numerical method to 
solve directly the elliptic problem Eq. I|2.5|l . and use the 
linearized solution as an analytic solution against which 
the numerical method can be validated. The reason for 
doing this is that we expect, in future work, to apply the 
numerical method to problems for which an approximate 
analytic solution cannot be found. 

A. Linearized Analysis 

If the perturbation introduced in Eq. I|4.1|) is small 
enough, i.e. e <C 1, it is a good approximation to neglect 
terms of order e 2 and higher. In Sec. IIVBI we consider 
how small e should be. 

It turns out that the ansatz 

$ = 1/r, V(r,6) = V s (r)+ew n (r)P n (cos6) (4.2) 

separates the linearized Hamiltonian constraint and we 
are left with an ordinary differential equation for w n (r) 

di<(r) + d 2 w' n {r) + d 3(n) w n (r) = (4.3) 

The coefficients were found using computer algebra 

di = r 2 {-l + 2V s {r)) 2 

d 2 =-r (_i + 2K s (r)) (3 - 6 V s (r) + 7r V&(r)) 
ds(n) =1 - ^H±ll + ( 3n (n + 1) - 8) V s (r) 
- 2 (3n(n + 1) - 10) V s (rf 
+ 4 (n(rc+l)-4) V s (rf + 7 r 2 Vg{rf . 

The fact, that the Schwarzschild solution Vs(r) 
(Eq. (|3.4[) ~) is itself a solution of the Hamiltonian con- 
straint, has been used. We will define the functions w n (r) 
to be those solutions of 1)4. 3|l which satisfy the boundary 
conditions 

w n (r = -) = !, Wn (r -> oo) = 0, (4.4) 



(although, in numerical calculations, it will be convenient 
to impose the outer boundary condition at finite r). 

Of course, a general solution to the perturbed 
Schwarzschild problem can be constructed by summing 
the eigenfunctions, i.e. 

oo 

V{r,9) = V s (r) + eJ2anW n (r)Pn{cos9) (4.5) 

n=l 

where the a n are arbitrary constants. 

1. The York tensor 

In order to show that the 3-metric defined by Eq l|4.2|) 
is not conformally flat, we need to prove that the York 
tensor fl4L fl5| 

Yijk = Rij-k - Rik;j + -l{R,j9ik — R,k9ij), (4.6) 

does not vanish identically. We have used computer alge- 
bra to find Yijk'. if only the first eigenfunction is present 
(i.e., if n = 1) then the York tensor is zero to first order 
in e; but, for example, 

II12 - e 2~n(n + l) w , 
4r 

So it is proven that the eigenfunctions with n > 1 do not 
represent a conformally flat geometry. 

B. Numerical Computation 

1. Implementation 

We want to solve Eq. 12.5(1 . subject to the boundary 
condition defined by Eq. (|4.1|l and 

V(r,6)\ r=roB =V 3 (r) 

with Vs(r) given by Eq. I|3.4|l . with a standard numerical 
elliptic solver. The computations were done using the 
Cactus-Computational- Toolkit 0] and the TAT-Jacobi 
elliptic solver 01 using the Jacobi method pj| imple- 
menting Eq. (|2.5|l and the boundary conditions by means 
of second order finite differencing. 

Due to the symmetry of the problem we impose the 
additional boundary conditions 



where n is normal to the 9 = 0, 7r-surfaces. For simplicity 
we implemented these conditions at 9 = + r;, 7r — r/ with 
r) <C 1 and of the order of magnitude of the accuracy up 
to which Eq. I|2.5|) is to be solved. This is a common 
procedure to avoid numerical problems at the singular 



points of the equation (9 = 0, w) - for example, the same 
situation arises in the case of the Laplace equation in 
spherical coordinates. 

To solve non-linear elliptic PDEs using the Jacobi 
method the choice of the initial guess for V(r, 9) is cru- 
cial. It turns out to be sufficient to use Eq. (|3.4H plus 
a Legendre perturbation of a given order n with linearly 
decreasing amplitude from the inner to the outer bound- 
ary consistent with the boundary conditions above. To 
compute Vs(r) for the initial guess numerically, Eq. (|3.3|) 
was solved for r by a numerical integration of Eq. (|3.2() 
which was then substituted into Eq. (|3.4|l - the problem 
being that Eq. 1)3.3)1 is explicit in the wrong direction. 

The convergence of TAT- Jacobi is very slow mainly 
due to non-linear terms in Eq. 1)2.5)1 which are dominant 
close to the horizon. So we first run the elliptic solver 
with the residual multiplied by (1 — 2V), which is small 
near the horizon. This gives the solver the opportunity 
to get an accurate solution everywhere else before, in a 
second run, we solve the original equation. By means 
of this technique the convergence speed was significantly 
increased. 



2. Results 

We investigate the case n — 2 with the outer boundary 
at roB — 10 and e = 0.005. We show that three dif- 
ferent resolutions exhibit second order convergence, and 
compare the numerical results to the linearized ones of 
Sec. IIV Al The next table describes the different resolu- 
tions used. With final residual we mean the residual of 
Eq. H2.5fl after the elliptic solver has finished. 



N r 


N B A r Ag final 


resid. 


of Eq. 


low 97 


33 0.098 0.095 


4.0 


10-'' 


medium 193 


65 0.049 0.048 


1.0 


io- 7 


high 385 


129 0.025 0.024 


2.5 


10" 8 



Fig. n shows the full numerical results together with 
the linearized one at 9 = 0.59. In the numerical case, 
we plot the difference V(r, 9) — Vs(r) and normalized to 
unity at r = \ . The graphs suggest second order conver- 
gence which is confirmed by Fig. El It is obvious that the 
linearized solution is not the limit of infinite resolution. 

We have thus shown that our numerics show the right 
convergence, but it must still be understood why the 
linearized solution deviates so much from the numeri- 
cal ones. The answer is that higher order terms get large 
close to the horizon. Using computer algebra, we substi- 
tuted the linearized solution which fulfills the linearized 
Hamiltonian constraint Eq. 1)4. 3J) up to an error of 10 -6 
into Eq. 1)2. 5JI . The residual is shown in Fig. and we 
also show 100 x the residual when e is smaller by a factor 
of 10, i.e. e = 0.0005. Although the linearized solution 
is not accurate close to the boundaries for e = 0.005, the 
residual shows the right quadratic scaling for decreasing 
e due to second order terms. This gives us a measure 
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FIG. 1: Linearized and full numerical solutions for n — 2, 
9 = 0.59 
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FIG. 2: Confirmation of second order convergence in the case 
n = 2, 9 = 0.59 

of how small e should be in order to obtain a desired 
accuracy. 

The computations were repeated for different positions 
of the boundaries and different n, recovering convergence 
as before. Fig. 21 shows the linearized solutions - i.e. full 
solutions provided e is small enough - for different n. 



V. CONCLUSIONS 

The paper has presented a new method of finding ini- 
tial data for the Einstein equations. The method is not 
based upon a conformal decomposition, but instead a 
unit vector field is chosen and then the Hamiltonian con- 
straint reduces to a partial differential equation in one 
unknown. At least in the case considered here - zero 
extrinsic curvature, axisymmetry, and with the vector 
field being a normalized gradient of a scalar field - the 
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FIG. 3: Residuum of Eq. I)2.5fl for the linearized solution 
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FIG. 4: Linearized solutions for different n 



At the analytic level, we have constructed a coordi- 
nate transformation from isotropic Schwarzschild coordi- 
nates into coordinates in which the metric satisfies our 
ansatz. Further, we have investigated small perturba- 
tions of Schwarzschild, and have shown that in this case 
the Hamiltonian constraint is separable so that its so- 
lution may be written as an infinite sum of eigenfunc- 
tions. The York tensor of the perturbed solution has 
been found: except for the case of the first eigenfunction, 
the York tensor is non-zero and thus the 3-metric is not 
conformally flat. 

Two immediate items of further work are envisaged 

• In the perturbed Schwarzschild, case it will be in- 
teresting to evolve the initial data obtained here 
and compare it to an evolution from initial data 
obtained from the standard conformal decomposi- 
tion. 

• The two black hole problem needs to be tackled. 
This will be difficult technically because at some 
point between the two black holes the unit vector 
field must be singular, which means that some of 
the coefficients in Eq. (12.51) will be singular. At the 
analytic level this simply means that V — at the 
singular point, but the existence of a singularity 
will complicate the numerics. 

The conformal method for finding initial data has been 
known and refined over many years. Our method re- 
quires further development - and in particular a solution 
to the two black hole problem - before some form of sys- 
tematic comparison between the two methods would be 
worthwhile. Even so, the existence of an alternative to 
the conformal method may be useful in investigating how 
different initial data sets, representing the same physics, 
affect spacetime evolution. 



Hamiltonian constraint is elliptic, and we have developed 
a numerical programme for solving this equation. Wc 
have tested our numerical routine with boundary data 
corresponding to both the unperturbed and perturbed 
Schwarzschild solutions, and the observed second order 
convergence is a validation of the numerical method. 
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